function [handle_vec,fmat,xmat,ymat]=plot_surfunc(funcmod,xgrid,ygrid,xlabel,ylabel,varargin)
% ========================================================================
% Given function handle *funcmod* and a grid of values for x and y, *xgrid*
% and *ygrid* evaluate the function at the mesh of the two grids and plot 
% the surface as well as the contours
% 
% This function can easily be adapted to DSGEs using the Kalman Filter 
% Alejandro Justiniano 
% ========================================================================

[xmat,ymat]=meshgrid(xgrid,ygrid); 
nx=length(xgrid); 
ny=length(ygrid); 
fmat=zeros(nx,ny); 
count=0;
for hh=1:nx;    
    for ii=1:ny;
        count = count + 1;
        fmat(hh,ii)=feval(funcmod,xgrid(hh),ygrid(ii),varargin{:});
    end
end
xmat=xmat'; 
ymat=ymat'; 

handle_vec=zeros(2,1); 

handle_vec(1)=figure; 
surf(xmat,ymat,fmat); 
YLabel(ylabel); 
XLabel(xlabel); 

handle_vec(2)=figure; 
CC=contour(xmat,ymat,fmat,15); 
YLabel(ylabel); 
XLabel(xlabel); 
clabel(CC); 